Este R Markdown recoge el enunciado de la práctica de la asignatura de redes sociales.

El objetivo es analizar un grafo, que se provee como fichero en el mismo paquete que este enunciado. En este fichero, encontramos solamente dos columnas, correspondiente a una interacción entre dos nodos de la red. Esta red está formada por distintos individuos que tienen contactos cara a cara durante un período de tiempo.

A continuación, dividimos la práctica en apartados, con una breve descripción de qué debe contener cada chunk de código donde el alumno desarrollará su respuesta así como las explicaciones que considere oportunas. Por favor, razona todas tus soluciones y escribe las explicaciones en azul.

Junto al título de cada apartado se encuentra la puntuación del mismo (pueden obtenerse hasta 10,5 puntos, aunque solamente se evaluará del 0 al 10).

Carga de datos y comprobaciones iniciales (0,5 puntos)

En este apartado, se pide:

Cargar el fichero adjunto en la práctica.

library(igraph)

rm(list = ls())

setwd('C:/Users/Javier/Documents/masterAFI/18. Analisis de Grafos/01. Redes sociales/practica/')
dd <- read.csv('red_contactos.csv', sep = ';')

Convertirlo en un objeto grafo de IGraph.

gg <- graph.data.frame(dd, directed = FALSE)
summary(gg)
## IGRAPH ea482fa UN-- 1390 222744 -- 
## + attr: name (v/c)

Comprobar que, efectivamente, tiene el número de nodos y enlaces correcto.

vcount(gg)
## [1] 1390
ecount(gg)
## [1] 222744

Simplificar: eliminar bucles y agregar enlaces múltiples

gg2 <- simplify(gg, remove.multiple = TRUE, remove.loops = TRUE)

# para cada enlace, calculamos cuantos caminos cortos hay entre sus vértices 
# el camino mínimo será igual al enlace, por lo tanto, la cantidad será igual al número de enlaces múltiples
# esta cantidad, lo añadimos como peso

E(gg2)$weight = sapply(E(gg2), function(e) {
  length(all_shortest_paths(gg, from = ends(gg2, e)[1], to = ends(gg2, e)[2])$res) } )


summary(gg2)
## IGRAPH ea5c390 UNW- 1390 53942 -- 
## + attr: name (v/c), weight (e/n)
hist(E(gg2)$weight, breaks = 30)

Selección de la componente conexa mayor (0,5 puntos)

En este apartado, se pide realizar los pasos adecuados para generar un nuevo objeto grafo, que sea conexo, y que involucre a todos los nodos y enlaces de la componente conexa mayor del grafo original.

Comprobamos que no sea conexo:

is.connected(gg2)
## [1] FALSE

Visualizamos los tamaños de las distintas componentes conexas:

ccs <- components(gg2)
ccs$csize
## [1] 1388    1    1

Vemos que hay 3 componentes conexas, una de ellas acumula prácticamente todos los puntos. Nos quedamos con ella:

id_compmayor <- which.max(ccs$csize)

vids <- ccs$membership == id_compmayor
vids <- which(ccs$membership == id_compmayor)

gg3 <- induced_subgraph(gg2, vids = vids )

summary(gg3)
## IGRAPH 1ae7aa5 UNW- 1388 53942 -- 
## + attr: name (v/c), weight (e/n)

Análisis descriptivo de la componente conexa mayor (2,5 puntos)

En este apartado, se pide analizar descriptivamente el grafo usando los conceptos que hemos visto durante las clases de teoría:

Grado medio

mean(degree(gg3))
## [1] 77.72622
hist(degree(gg3), breaks = 20)

Distancia media

average.path.length(gg3)
## [1] 3.066029

Diámetro

diameter(gg3)
## [1] 14

Distribución de grados y ajuste a una Power-Law

Podemos empezar por ver la densidad:

graph.density(gg3)
## [1] 0.0560391

A continuación, ajustamos una Power Law:

deg_dist <- degree_distribution(gg3)
fit <- fit_power_law(deg_dist)
fit
## $continuous
## [1] TRUE
## 
## $alpha
## [1] 3.447821
## 
## $xmin
## [1] 0.005043228
## 
## $logLik
## [1] 391.6632
## 
## $KS.stat
## [1] 0.1568648
## 
## $KS.p
## [1] 0.03535439
plot(deg_dist + 0.0001, log = 'xy', xlab = 'Node Degree', ylab = 'Probability')
lines(seq(deg_dist), seq(deg_dist)^-fit$alpha, col='red')

Tanto con el p-value como con el plot, podemos ver que nuestros datos no se ajustan a una Power Law.

Podemos tratar de ver que pasaría si eliminamos aquellos nodos con los que obteníamos un valor 0 en la función de distribución:

gg_prueba <- induced_subgraph(gg3, vids = V(gg3)[deg_dist > 0])
deg_dist_prueba <- degree_distribution(gg_prueba)
fit_prueba <- fit_power_law(deg_dist_prueba)
fit_prueba
## $continuous
## [1] TRUE
## 
## $alpha
## [1] 3.387724
## 
## $xmin
## [1] 0.008536585
## 
## $logLik
## [1] 185.4567
## 
## $KS.stat
## [1] 0.1363636
## 
## $KS.p
## [1] 0.386501
plot(deg_dist_prueba  + 0.0001, log = 'xy', xlab = 'Node Degree', ylab = 'Probability')
lines(seq(deg_dist_prueba), seq(deg_dist_prueba)^-fit_prueba$alpha, col='red')

Con el nuevo grafo que hemos creado parece que sí que tendríamos unos datos que se pueden ajustar a una Power Lab.

Clustering

Podemos observar el número de triángulos por nodo:

triangulos <- count_triangles(gg3)
mean(triangulos)
## [1] 1037.44
hist(triangulos, breaks = 30)

Entropía de Shannon de los nodos

hist(diversity(gg3), breaks = 20, main = '', xlab = 'Diversity')

Centralidad de los nodos y comparación con métricas de grado y clustering

Vemos las distintas centralidades, hacemos histogramas de ellas y observamos como se correlacionan:

Degree <- degree(gg3)
Eig <- evcent(gg3)$vector
Closeness <- closeness(gg3)
Betweenness <- betweenness(gg3)

par(mfrow=c(2,2))
for (centrality in c('Degree', 'Eig', 'Closeness', 'Betweenness')){
  hist(get(centrality), breaks = 30, xlab = centrality, main = '')
}

centralities <- cbind(Degree, Eig, Closeness, Betweenness, triangulos)
round(cor(centralities), 2)
##             Degree  Eig Closeness Betweenness triangulos
## Degree        1.00 0.93      0.79        0.85       0.92
## Eig           0.93 1.00      0.69        0.78       0.92
## Closeness     0.79 0.69      1.00        0.59       0.59
## Betweenness   0.85 0.78      0.59        1.00       0.87
## triangulos    0.92 0.92      0.59        0.87       1.00

Podemos observar como las centralidades por grado, cercanía y betweeness se parecen bastante entre ellas, además de con el número de triángulos de cada nodo. Sin embargo, parece que Closeness es la que más difiere de las demás.

Análisis de comunidades de la componente conexa mayor (1,5 puntos)

En este apartado, se pide aplicar dos algoritmos de detección de comunidades, compararlos y seleccionar cuál es, en tu opinión, el que da una mejor respuesta. Razona tu selección.

Usaremos los algoritmos de Infomap y de Walktrap para detectar comunidades, calculamos las modularidades para cada uno y comparamos los resultados de los algoritmos:

comms_infomap <- infomap.community(gg3)
comms_walktrap <- walktrap.community(gg3)

modularity(comms_infomap)
## [1] 0.4686043
modularity(comms_walktrap)
## [1] 0.4634219
compare(comms_infomap, comms_walktrap, method='nmi')
## [1] 0.8131684

Vemos que la función compare nos devuelve un valor de 0.81, es decir, los algoritmos nos devuelven resultados más o menos parecidos. Para elegir uno de ellos, veremos como éstos dividen a sus comunidades:

table(comms_infomap$membership)
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   9  55 392  45  26  11  37  25  12  30  61  87  22  52  61  60  35   4  13  18 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##  31  35   1  17  38  13   3  25   6  41  10  22  14   3  30  13   7   8   6   6 
##  41  42 
##   2   2
table(comms_walktrap$membership)
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##  14  12  42  13 134 396  29  50  53  54  28  46   8  15  17  97   2   4  52  83 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##   5  10  16   7  30  13  19  32  17   2   2   1   1   1   1   1   1   1   1   1 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1

Vemos que el algoritmo de Infomap divide comunidades de una forma más uniforme, mientras que, el algoritmo Walktrap nos devuelve bastantes comunidades formadas por un solo miembro. Con estos resultados, decidimos quedarnos con Infomap como algoritmo de detección de comunidades.

Visualización del grafo por comunidades de la componente conexa mayor (1,5 puntos)

En este apartado, se pide visualizar el grafo coloreando cada nodo en función de la comunidad a la que pertenezca, según tu elección del apartado anterior.

wws <- ifelse(crossing(comms_infomap, gg3), 1, 500)
ll <- layout_with_fr(gg3, weights=wws)

palcol <- rainbow(n=max(comms_infomap$membership))

plot(gg3,
     layout = ll,
     vertex.label='',
     vertex.size=log(degree(gg3)),
     vertex.color=palcol[comms_infomap$membership])

Difundiendo un rumor (o un virus) en la componente conexa mayor (4 puntos)

Este apartado es el que más peso en la práctica tiene. Vamos a implementar un modelo epidemiológico sobre el grafo que, típicamente, se utiliza para simular escenarios de difusión de enfermedades pero también en contextos como la distribución de rumores e información. Vamos a implementar un modelo SIR que se caracteriza por tener los siguientes parámetros:

Se pide desarrollar una función que tenga como parámetros los tres valores anteriores y un cuarto que sea un grafo que, en nuestro caso, será la componente conexa mayor del grafo original de esta práctica. Dicha función simulará el proceso SIR:

Se pide ejecutar una simulación para tres o cuatro valores del parámetro beta (N y gamma pueden ser fijos en estas simulaciones) de este proceso de manera que se pueda visualizar:

En este primer bloque de código, vamos a definir las siguientes funciones:

generate_sir <- function(N = 5, beta = 0.1, gamma = 0.1, graph){
  graph <- set_vertex_attr(graph, 'label', value = 'S')
  graph <- set_edge_attr(graph, 'label', value = '')
  
  sample_index <- sample(1:vcount(graph), N)
  graph <- set_vertex_attr(graph, 'label', index = sample_index, value = 'I')
  
  n_contagiados = N
  n_susceptibles = vcount(graph) - N
  n_recuperados = 0
  iteraciones_sin_contagios = 0
  iteracion= 1
  
  lista_contagiados = c(n_contagiados)
  lista_susceptibles = c(n_susceptibles)
  lista_recuperados = c(n_recuperados)
  lista_nuevos_contagios = c(n_contagiados)
  lista_grafo = list(graph)
  
  while (iteraciones_sin_contagios < 3){
    n_nuevos_contagios = 0
    initial_graph <- graph
    
    for (v in V(initial_graph)){
      state <- V(initial_graph)[v]$label
      
      if (state == 'I'){
        if (rbinom(n=1, size=1, prob = gamma) == 1){
          V(graph)[v]$label <- 'R'
          n_recuperados = n_recuperados + 1
          n_contagiados = n_contagiados - 1
        }
      }
      if (state == 'S'){
        for (neighbor in neighbors(initial_graph, v)){
          if ((V(initial_graph)[neighbor]$label == 'I') && (rbinom(n=1, size=1, prob = beta) == 1)){
            V(graph)[v]$label <- 'I'
            E(graph)[v %--% neighbor]$label <- 'I'
            n_contagiados = n_contagiados + 1
            n_nuevos_contagios = n_nuevos_contagios + 1
            n_susceptibles = n_susceptibles - 1
            break
          }
        }
      }
    }
    if (n_nuevos_contagios == 0){
      iteraciones_sin_contagios = iteraciones_sin_contagios + 1
    } else {
      iteraciones_sin_contagios = 0
    }
    iteracion = iteracion + 1
    lista_contagiados = c(lista_contagiados, n_contagiados)
    lista_susceptibles = c(lista_susceptibles, n_susceptibles)
    lista_recuperados = c(lista_recuperados, n_recuperados)
    lista_nuevos_contagios = c(lista_nuevos_contagios, n_nuevos_contagios)
    lista_grafo[[iteracion]] = graph
  }
  
  return (list('contagiados' = lista_contagiados,
               'susceptibles' = lista_susceptibles,
               'recuperados' = lista_recuperados,
               'nuevos_contagios' = lista_nuevos_contagios,
               'iteraciones' = iteracion,
               'grafos' = lista_grafo
          ))

}

plot_evolution <- function(evolution){
  plot(1:evolution$iteraciones, evolution$contagiados, 
       type = 'l', col = 'red',
       ylab = 'Nº Nodos', xlab = 'Iteraciones', 
       ylim = c(0, max(evolution$susceptibles)),
       main = 'Evolución')
  lines(1:evolution$iteraciones, evolution$recuperados, 
        type = 'l', col='green')
  lines(1:evolution$iteraciones, evolution$susceptibles, 
        type = 'l', col='blue')
  legend(evolution$iteraciones, max(evolution$susceptibles)/2, 
         legend=c('Contagiados', 'Recuperados', 'Susceptibles'),
         col=c('red', 'green', 'blue'), 
         lty=1, cex=0.8, xjust = 1)
}

plot_new_infections <- function(list_new_infections, iteraciones){
  plot(1:iteraciones, list_new_infections, 
       type = 'o', col = 'red', log = 'y', 
       ylab = 'Nº Nuevos Infectados', xlab = 'Iteraciones',
       main = 'Nuevos infectados')
}

plot_infection_graph <- function(graph, index){
  graph_colors <- ifelse(V(graph)$label == 'I', 'red',
                       ifelse(V(graph)$label == 'R', 'green', 'blue'))

  subg <- subgraph.edges(graph, E(graph)[label == 'I'])
  
  plot(subg,
       layout = ll,
       vertex.label = '',
       vertex.size = log(degree(graph)),
       edge.label = '',
       vertex.color = graph_colors,
       main = paste('Iteración', index))
  
}

create_gif <- function(list_graphs){
  img <- magick::image_graph(width = 1000, height = 1000)
  for (i in 1:length(list_graphs)){
    plot_infection_graph(list_graphs[[i]], i)
  }
  dev.off()
  animation <- magick::image_animate(img, fps = 2, optimize = TRUE)
  print(animation)
}

Veremos distintas ejecuciones del simulador del proceso SIR, dibujando las gráficas y los GIF correspondientes a cada una de ellas.

Nuestras distintas ejecuciones tendrán como parámetros los mismos valores para gamma (probabilidad de recuperación; gamma = 0.1), N (número de nodos iniciales infectados; N = 5) y el grafo (utilizaremos el grafo resultante de los apartados anteriores). El parámetro que vamos a ir variando será el beta.

Beta = 0.1

ejemplo_SIR_1 <- generate_sir(graph = gg3, beta = 0.1)

plot_evolution(ejemplo_SIR_1)

plot_new_infections(ejemplo_SIR_1$nuevos_contagios, ejemplo_SIR_1$iteraciones)

create_gif(ejemplo_SIR_1$grafos)
##    format width height colorspace matte filesize density
## 1     gif  1000   1000       sRGB  TRUE        0   72x72
## 2     gif   784    859       sRGB  TRUE        0   72x72
## 3     gif   442    687       sRGB  TRUE        0   72x72
## 4     gif   605    809       sRGB  TRUE        0   72x72
## 5     gif   620    857       sRGB  TRUE        0   72x72
## 6     gif   616    792       sRGB  TRUE        0   72x72
## 7     gif   620    857       sRGB  TRUE        0   72x72
## 8     gif   613    857       sRGB  TRUE        0   72x72
## 9     gif   627    857       sRGB  TRUE        0   72x72
## 10    gif   592    857       sRGB  TRUE        0   72x72
## 11    gif   590    857       sRGB  TRUE        0   72x72
## 12    gif   535    856       sRGB  TRUE        0   72x72
## 13    gif   587    859       sRGB  TRUE        0   72x72
## 14    gif   546    857       sRGB  TRUE        0   72x72
## 15    gif   769    856       sRGB  TRUE        0   72x72
## 16    gif   595    687       sRGB  TRUE        0   72x72
## 17    gif   754    695       sRGB  TRUE        0   72x72
## 18    gif   769    856       sRGB  TRUE        0   72x72
## 19    gif   244    682       sRGB  TRUE        0   72x72
## 20    gif   266    738       sRGB  TRUE        0   72x72
## 21    gif   490    653       sRGB  TRUE        0   72x72
## 22    gif   349    688       sRGB  TRUE        0   72x72
## 23    gif   769    857       sRGB  TRUE        0   72x72
## 24    gif   241    687       sRGB  TRUE        0   72x72
## 25    gif   750    781       sRGB  TRUE        0   72x72
## 26    gif   206    630       sRGB  TRUE        0   72x72
## 27    gif   426    578       sRGB  TRUE        0   72x72
## 28    gif   205    685       sRGB  TRUE        0   72x72

Beta = 0.05

ejemplo_SIR_2 <- generate_sir(graph = gg3, beta = 0.05)

plot_evolution(ejemplo_SIR_2)

plot_new_infections(ejemplo_SIR_2$nuevos_contagios, ejemplo_SIR_2$iteraciones)

create_gif(ejemplo_SIR_2$grafos)
##    format width height colorspace matte filesize density
## 1     gif  1000   1000       sRGB  TRUE        0   72x72
## 2     gif   784    859       sRGB  TRUE        0   72x72
## 3     gif   443    786       sRGB  TRUE        0   72x72
## 4     gif   612    809       sRGB  TRUE        0   72x72
## 5     gif   626    792       sRGB  TRUE        0   72x72
## 6     gif   613    792       sRGB  TRUE        0   72x72
## 7     gif   619    809       sRGB  TRUE        0   72x72
## 8     gif   627    792       sRGB  TRUE        0   72x72
## 9     gif   614    856       sRGB  TRUE        0   72x72
## 10    gif   620    856       sRGB  TRUE        0   72x72
## 11    gif   614    857       sRGB  TRUE        0   72x72
## 12    gif   613    856       sRGB  TRUE        0   72x72
## 13    gif   615    856       sRGB  TRUE        0   72x72
## 14    gif   599    856       sRGB  TRUE        0   72x72
## 15    gif   592    857       sRGB  TRUE        0   72x72
## 16    gif   592    856       sRGB  TRUE        0   72x72
## 17    gif   562    695       sRGB  TRUE        0   72x72
## 18    gif   554    856       sRGB  TRUE        0   72x72
## 19    gif   574    702       sRGB  TRUE        0   72x72
## 20    gif   592    857       sRGB  TRUE        0   72x72
## 21    gif   592    856       sRGB  TRUE        0   72x72
## 22    gif   594    856       sRGB  TRUE        0   72x72
## 23    gif   592    856       sRGB  TRUE        0   72x72
## 24    gif   286    681       sRGB  TRUE        0   72x72
## 25    gif   343    639       sRGB  TRUE        0   72x72
## 26    gif   592    856       sRGB  TRUE        0   72x72
## 27    gif   614    857       sRGB  TRUE        0   72x72
## 28    gif   524    854       sRGB  TRUE        0   72x72
## 29    gif   615    857       sRGB  TRUE        0   72x72
## 30    gif   467    694       sRGB  TRUE        0   72x72
## 31    gif   578    738       sRGB  TRUE        0   72x72
## 32    gif   361    692       sRGB  TRUE        0   72x72

Beta = 0.01

ejemplo_SIR_3 <- generate_sir(graph = gg3, beta = 0.01)

plot_evolution(ejemplo_SIR_3)

plot_new_infections(ejemplo_SIR_3$nuevos_contagios, ejemplo_SIR_3$iteraciones)

create_gif(ejemplo_SIR_3$grafos)
##    format width height colorspace matte filesize density
## 1     gif  1000   1000       sRGB  TRUE        0   72x72
## 2     gif   784    859       sRGB  TRUE        0   72x72
## 3     gif   441    628       sRGB  TRUE        0   72x72
## 4     gif   442    676       sRGB  TRUE        0   72x72
## 5     gif   473    693       sRGB  TRUE        0   72x72
## 6     gif   565    691       sRGB  TRUE        0   72x72
## 7     gif   545    693       sRGB  TRUE        0   72x72
## 8     gif   606    809       sRGB  TRUE        0   72x72
## 9     gif   600    792       sRGB  TRUE        0   72x72
## 10    gif   619    792       sRGB  TRUE        0   72x72
## 11    gif   620    792       sRGB  TRUE        0   72x72
## 12    gif   617    792       sRGB  TRUE        0   72x72
## 13    gif   620    792       sRGB  TRUE        0   72x72
## 14    gif   603    793       sRGB  TRUE        0   72x72
## 15    gif   613    792       sRGB  TRUE        0   72x72
## 16    gif   612    792       sRGB  TRUE        0   72x72
## 17    gif   620    792       sRGB  TRUE        0   72x72
## 18    gif   612    792       sRGB  TRUE        0   72x72
## 19    gif   613    792       sRGB  TRUE        0   72x72
## 20    gif   604    792       sRGB  TRUE        0   72x72
## 21    gif   587    809       sRGB  TRUE        0   72x72
## 22    gif   603    792       sRGB  TRUE        0   72x72
## 23    gif   620    792       sRGB  TRUE        0   72x72
## 24    gif   539    792       sRGB  TRUE        0   72x72
## 25    gif   603    792       sRGB  TRUE        0   72x72
## 26    gif   603    792       sRGB  TRUE        0   72x72
## 27    gif   568    792       sRGB  TRUE        0   72x72
## 28    gif   603    792       sRGB  TRUE        0   72x72
## 29    gif   520    710       sRGB  TRUE        0   72x72
## 30    gif   603    792       sRGB  TRUE        0   72x72
## 31    gif   532    703       sRGB  TRUE        0   72x72
## 32    gif   560    703       sRGB  TRUE        0   72x72
## 33    gif   604    792       sRGB  TRUE        0   72x72
## 34    gif   560    785       sRGB  TRUE        0   72x72
## 35    gif   435    652       sRGB  TRUE        0   72x72
## 36    gif   512    681       sRGB  TRUE        0   72x72
## 37    gif   393    682       sRGB  TRUE        0   72x72
## 38    gif   271    602       sRGB  TRUE        0   72x72
## 39    gif   384    558       sRGB  TRUE        0   72x72